A study of H+H2 and several H-bonded molecules by phaseless auxiliary-field 
quantum Monte Carlo with planewave and Gaussian basis sets 
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We present phaseless auxiliary-field (AF) quantum Monte Carlo (QMC) calculations of the ground 
states of some hydrogen-bonded systems. These systems were selected to test and benchmark dif- 
ferent aspects of the new phaseless AF QMC method. They include the transition state of H-I-H2 
near the equilibrium geometry and in the van der Walls limit, as well as H2O, OH, and H2O2 
molecules. Most of these systems present significant challenges for traditional independent-particle 
electronic structure approaches, and many also have exact results available. The phaseless AF QMC 
method is used either with a planewave basis with pseudopotentials or with all-electron Gaussian 
basis sets. For some systems, calculations are done with both to compare and characterize the 
performance of AF QMC under different basis sets and different Hubbard-Stratonovich decomposi- 
tions. Excellent results are obtained using as input single Slater determinant wave functions taken 
from independent-particle calculations. Comparisons of the Gaussian based AF QMC results with 
exact full configuration show that the errors from controlling the phase problem with the phaseless 
approximation are small. At the large basis-size limit, the AF QMC results using both types of 
basis sets are in good agreement with each other and with experimental values. 

PACS numbers: 



I. INTRODUCTION 



Quantum Monte Carlo (QMC) methods [1, 2] offer a 
unique way to treat explicitly the many-electron prob- 
lem. The many-body solution is obtained in a statisti- 
cal sense by building stochastic ensembles that sample 
the wave function in some representation. This leads to 
computational costs that scale as a low power with the 
number of particles and basis size. Although in prac- 
tice QMC methods are often not exact, they have shown 
considerably greater accuracy than traditional electronic 
structure approaches in a variety of systems. They are 
increasingly applied and are establishing themselves as 
a unique approach for studying both realistic materials 
and important model systems. 

Recently, a new phaseless auxiliary-field QMC (AF 
QMC) method has been developed and applied for elec- 
tronic structure calculations [2, 3]. This method is formu- 
lated in a many-particle Hilbert space whose span is de- 
fined by a single-particle basis set. The freedom to choose 
the basis set can potentially result in increased efficiency. 
This can be very useful both for quantum chemistry ap- 
plications and in calculations with model Hamiltonians. 
Further, it is straightforward in this method to exploit 
well-established techniques of independent-particle the- 
ories for the chosen basis set. The ability to use any 
single-particle basis is thus an attractive feature of the 
AF QMC method. On the other hand, the use of finite 
basis sets requires monitoring the convergence of calcu- 
lated properties and extrapolation of the results to the 
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infinite basis size limit. 

Planewave and Gaussian basis sets are the most widely 
used in electronic structure calculations. Planewaves are 
appealing, because they form a complete orthonormal 
basis set, and convergence with respect to basis size is 
easily controlled. A single energy cutoff parameter Ecnt 
controls the basis size by including all planewaves with 
wavevector k such that k'^/2 < Ecut (Hartree atomic 
units are used throughout the paper). The infinite basis 
limit is approached by simply increasing -Ecut [4] . Local- 
ized basis sets, by contrast, offer a compact and efficient 
representation of the system's wavefunctions. Moreover, 
the resulting sparsity of the Hamiltonians can be very 
useful in 0{M) methods. Achieving basis set conver- 
gence, however, requires more care. For Gaussian basis 
sets, quantum chemists have compiled lists of basis sets 
of increasing quality for most of the elements [5] . Some 
of these basis sets have been designed for basis extrapola- 
tion not only in mean-field theories, but also in correlated 
calculations [6, 7]. 

The AF QMC method also provides a different route to 
controlling the Fermion sign problem [2, 8-10]. The stan- 
dard difllusion Monte Carlo (DMC) method [1, 11, 12] 
employs the fixed-node approximation [11] in real coor- 
dinate space. The AF QMC method uses random walks 
in a manifold of Slater determinants (in which antisym- 
metry is automatically imposed on each random walker). 
The Fermion sign/phase problem is controlled approxi- 
mately according to the overlap of each random walker 
(Slater determinant) with a trial wave function. Applica- 
tions of the phaseless AF QMC method to date, includ- 
ing second-row systems [2] and transition metal molecules 
[13] with planewave basis sets, and first-row [3] and post- 
d [14] molecular systems with Gaussian basis sets, indi- 
cate that this often reduces the reliance of the results on 
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the quality of the trial wave function. For example, with 
single determinant trial wave functions, the calculated to- 
tal energies at equilibrium geometries in molecules show 
typical systematic errors of no more than a few milli- 
Hartrees compared to exact/experimental results. This is 
roughly comparable to that of CCSD(T) (coupled-cluster 
with single and double excitations plus an approximate 
treatment of triple excitations). For stretched bonds in 
H2O [3] as well as N2 and F2 [15], the AF QMC method 
exhibits better overall accuracy and a more uniform be- 
havior than CCSD(T) in mapping the potential energy 
curve. 

The key features of the AF QMC method are thus 
its freedom of basis choice and control of the Fermion 
sign/phase problem via a constraint in Slater determi- 
nant space. The motivation for this study is therefore 
two- fold. First, we would like to further benchmark the 
planewave AF QMC method in challenging conditions, 
with large basis sets and correspondingly many auxiliary 
fields. Here we examine the transition state of the H2+H 
system as well as several hydrogen-bonded molecules. 
These are relatively simple systems, which have been dif- 
ficult for standard independent-electron methods and for 
which various results are available for comparison. Sec- 
ondly, we are interested in comparing the performance of 
the AF QMC method using two very different basis sets, 
namely, planewave basis sets together with pseudopo- 
tentials, and all-electron Gaussian basis sets. For this, 
calculations are carried out with Gaussian basis sets for 
H2 and the H2-I-H transition state, and comparisons are 
made with the planewave calculations. Additional Gaus- 
sian benchmark calculations are carried out in collinear 
H2-I-H near the Van der Waals minimum, which requires 
resolution of the energy on extremely small scales. 

The rest of the paper is organized as follows. In the 
next section we outline the relevant formalism of the AF 
QMC method. The planewave and pseudopotential re- 
sults are presented in Sec. Ill, including a study of the 
dissociation energy and the potential energy curve of H2, 
the transition state of II3, and the dissociation energies 
of several hydrogen-bonded molecules. In Sec. IV, we use 
a Gaussian basis to study the potential energy curves of 
H2 and H-I-H2, and compare some of these results with 
the AF QMC planewave results. Finally, we conclude in 
Sec. V with a brief summary. 



II. AF QMC METHOD 

The auxiliary-field quantum Monte Carlo method has 
been described elsewhere [2, 3]. Here we outline the rel- 
evant formulas to facilitate the ensuing discussion. The 
method shares with other QMC methods its use of the 
imaginary-time propagator e~^^ to obtain the ground 
state |*g) of H: 

|*g) oc lim e-l^"\^T)- (1) 

/3— >oo 



The ground state is obtained by filtering out the excited 
state contributions in the trial wave function |\I't), pro- 
vided that I^'t) has a non-zero overlap with \^g)- 

The many-body electronic Hamiltonian H can be writ- 
ten in any one-particle basis as, 

H = ffi +F2; 

-^2 = ^ ^ Vijklc\a<^]^^,Ck,a'Cl,a, (2) 

where and Ci,o- are the corresponding creation and 
annihilation operators of an electron with spin a in the 
i-th orbital (size of single-particle basis is M). The one- 
electron and two-electron matrix elements (Ty and Vijki ) 
depend on the chosen basis, and are assumed to be spin- 
independent. 

Equation. (1) is realized iteratively with a small time- 
step r such that (3 = Nt, and the /3 ^ oo limit is 
realized by letting N ^ oo. In this case, the Trot- 
ter decomposition of the propagator e~'^^: e~'^^ = 
g-r_f/i/2g-r_f/2p-Tffi/2 _^(^^^3^ Icads to Trottcr time-step 
errors, which can be removed by extrapolation, using sep- 
arate calculations with different values t. 

The central idea in the AF QMC method is the use of 
the Hubbard-Stratonovich (HS) transformation [16]: 

=n ^«"). (3) 

to map the many-body problem exemplified in H2 onto a 
linear combination of single-particle problems using only 
one-body operators The full many-body interaction 
is recovered exactly through the interaction between the 
one-body operators {va}, and all of the external auxiliary 
fields {aa}- This map relies on writing the two-body 
operator in a quadratic form, such as 

^2 = -^Ec«*- (4) 

a 

with a real number. This can always be done, as we 
illustrate below using first a planewave basis, and then 
any basis set. 

In a planewave basis set, the electron-electron interac- 
tion operator H2 can be written as; 

^2 = ^2^ Z^-r2-Ck+q,aCk'-q,a'Ck',a'Ck,a 
k,k',(T,a-' q^O ^ 

= i E ^ [/3(q)p(-q) + h.c.] + H[. (5) 

q>0 ^ 

Here cj^ ^ and Ck^o- are the creation and annihilation op- 
erators of an electron with momentum k and spin a. Q is 
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the supercch volume, k and k' are plancwaves within the 
cutoff radius, and the q vectors satisfy |k-|-qp/2 < E^ut- 
'^k+q.o-'^k.cr is a Fourier component of the 
electron density operator, and H[ is a one-body term 
which arises from the reordering of the creation and anni- 
hilation operators. For each wavevector q, the two-body 
term in the final expression in Eq. (5) can be expressed in 
terms of squares of the one-body operators proportional 
to /5(q) -I- /5(— q) and /5(q) — q), which become the 
one-body operators iice in Eq. (4). 

An explicit HS transformation can be given for any 
general basis as follows (more; efficient transformations 
may exist, however). The two-body interaction ma- 
trix Vijki is first expressed as a Hermitian supermatrix 
Vf^[i^i],iy[k,j] where yL,v = l,...,Af^. This is then ex- 
pressed in terms of its eigenvalues (—Ac) and eigenvectors 
V^,- = -Ea^«^*"^". The two-body operator 
H'2 of Eq. (2) can be written as the sum of a one-body 
operator H[ and a two-body operator H!^. The latter can 
be further expressed in terms of the eigenvectors oiV^^i, 
as 
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4 ^ 



Aq. 



where the one-body operators A^ are defined as 



(6) 



(7) 



Similar to the planewave basis, for each non-zero eigen- 
value Aq,, there are two one-body operators propor- 
tional to Aq, + AJj and Aq — Aj^. If the chosen basis set is 
real, then the HS transformation can be further simpli- 
fied, and the number of auxiliary fields will be equal to 
only the number of non-zero eigenvalues Aq [3] . 

The phasclcss AF QMC method [2] used in this paper 
controls the phase/sign problem [2, 10] in an approximate 
manner. The method recasts the imaginary-time path in- 
tegral as branching random walks in Slater-determinant 
space [10]. It uses a trial wave function to construct 
a complex importance-sampling transformation and to 
constrain the paths of the random walks. The ground- 
state energy, computed with the so-called mixed estima- 
tor, is approximate and not variational in the phaseless 
method. The error depends on \^t), vanishing when 
I^'t) is exact. This is the only imcontroUcd error in the 
method, in that it cannot be eliminated systematically. 
In applications to date, |^'t) has been taken as a sin- 
gle Slater determinant directly from mean-field calcula- 
tions, and the systematic error is shown to be quite small 
[2, 3, 13, 14]. 



III. RESULTS USING PLANEWAVE BASIS 

SETS 

Plancwaves arc more suited to periodic systems and 
require pseudopotentials to yield a tractable number of 



TABLE I: Planewave based calculations of the binding en- 
ergy of H2 vs. supercell size. DFT/GGA and the phaseless 
AF QMC results are shown. All energies are in eV, and su- 
percell dimensions are in atomic units. For comparison, the 
all-electron GGA number is 4.568 eV [17]. Statistical errors 
are on the last digit and are shown in parenthesis. The exact 
theoretical value is 4.746 eV [18], and the experimental value 
is 4.75 eV (with zero-point energy removed). 



supcrct'll 



dft/(;(;a 



AF QMC 



11x9x7 

12x10x9 

14x12x11 

16x12x11 

22x18x14 

00 



4.283 
4.444 
4.511 
4.512 
4.530 
4.531 



4.36(1) 
4.57(1) 
4.69(1) 
4.70(1) 
4.74(2) 



basis functions. However, isolatc;d niolec;ules can be stud- 
ied with plancwaves by employing periodic boundary 
conditions and large supercells, as in standard density 
functional theory (DFT) calculations. This is disad- 
vantageous, because one has to ensure that the super- 
cells are large enough to control the spurious interac- 
tions between the periodic images of the molecule. For 
a given planewave cutoff energy i?cut, the size of the 
planewave basis increases in proportion to the volume 
of the supercell. Consequently, the computational cost 
for the isolated molecule tends to be higher than using 
a localized basis, as we hirther disc;uss in Sec. IV. Al- 
though the planewave basis calculations are expensive, 
they are nevertheless valuable as they show the robust- 
ness and accuracy of the phasclcss AF QMC method 
for extremely large basis sets (and correspondingly many 
auxiliary fields). 

Here we study H2, H3, and several other hydrogen- 
bonded molecules H2O, OH, and H2 02- As is well known, 
first-row atoms like oxygen are challenging, since they 
have strong or "hard" pseudopotentials and require rela- 
tively large planewave basis sets to achieve convergence. 
Even in hydrogen, where there are no core electron states, 
pseudopotentials are usually used, since they significantly 
reduce the planewave basis size compared to treating 
the bare Coulomb potential of the proton. The hydro- 
gen and oxygen pseudopotentials are generated by the 
OPIUM program [19], using the neutral atoms as refer- 
ence configurations. The cutoff radii used in the gen- 
eration of the oxygen pseudopotentials arc rc{s) = 1.05 
and rc{p) = 1.02 Bohr, where s and p correspond to 
I = and I = 1 partial waves, respectively. For hydrogen 
i'c{s) = 0.66 Bohr was used. These relatively small re's 
are needed for both atoms, due to the short bondlengths 
in H2O, and result in relatively hard pseudopotentials. 
Small re's, however, generally result in pseudopotentials 
with better transferability. In all of the studies shown 
below, the same pseudopotentials were used, even in 
molecules with larger bondlengths. The Ecut needed with 
these pseudopotentials is about 41 Hartrce. This Ecut 
was chosen such that the resulting planewave basis con- 



4 




R (a.u.) 



FIG. 1: The potential energy curve of H2 as obtained by AF 
QMC with a planewave basis and a hydrogen pseudopotential. 
We show also a Morse potential fit for the QMC data. The 
QMC equilibrium bondlength from the fit is 1.416(4) Bohr to 
be compared with the exact value is 1.40083 Bohr. Supercell 
used is 16 X 12 X 11 Bohr^. 



vergence errors are less than a few meV in DFT calcu- 
lations. A roughly similar planewave basis convergence 
error is expected at the AF QMC level, based on previous 
applications in TiO and other systems [13, 20, 21]. These 
convergence errors are much smaller than the QMC sta- 
tistical error. 

The quality of the pseudopotentials is further assessed 
by comparing the pseudopotential calculations with all- 
electron (AE) results using density functional methods, 
which tests the pseudopotentials, at least at the mean- 
field level. In all of the cases reported in this study, we 
found excellent agreement between AE and pseudopoten- 
tial results, except in cases where the non-linear core cor- 
rection error [22, 23] is important in DFT (all molecules 
containing oxygen), as we discuss in Sec. IIIC. 

As mentioned, AF QMC relies on a trial wave func- 
tion to control the phase problem. In the planewave cal- 
culations, we used a single Slater determinant from a 
planewave based density functional calculation obtained 
with a GGA functional [24], with no further optimiza- 
tions. 



A. Ortho- and Para-H2 molecule 

Table I summarizes results for the binding energy of 
the H2 molecule, using DFT/GGA and AF QMC for 
several supercells, and compares these to exact results 
[18] and experiment. The experimental bondlength of 
H2 was used in all the calculations. The binding en- 
ergy is calculated as the difference in energy between the 
H atom (times two) and the molecule, each placed in 
the same supercell. The density functional binding en- 
ergy obtained using the hydrogen pseudopotential con- 



TABLE II: Planewave based AF QMC energies of the ortho- 
(^E) and para-H2 ("^E) molecule for two supercell sizes. The 
bondlength was fixed at ii = 1.42 Bohr in all cases. All ener- 
gies are in eV. The exact calculated energy gap A is 10.495 eV 
[18]. Statistical errors are on the last digit and are shown in 
parenthesis. 



supercell 




3E 


A 


11x9x7 


-32.59(1) 


-22.329(4) 


10.26(1) 


22x18x14 


-32.01(2) 


-21.546(7) 


10.46(2) 



verges, with respect to size effect, to 4.531 eV, which is 
in reasonable agreement with the all-electron (i.e. us- 
ing the proton's bare Coulomb potential) binding energy 
4.568 eV obtained using NWCHEM [17], and with the 
all-electron value of 4.540 eV reported in Ref. [25]. The 
agreement between the pseudopotential and all-electron 
results is a reflection of the good transferability of the 
hydrogen pseudopotential. The AF QMC binding en- 
ergy with the largest supercell is 4.74(2) eV, which is 
in excellent agreement with the experimental value of 
4.75 eV (zero point energy removed) and the exact cal- 
culated value of 4.746 eV [18]. 

Figure 1 shows the H2 AF QMC potential energy 
curve, using a 16 x 12 x 11 Bohr'^ supercell. Finite size 
effects, as in Table I, likely vary with the H2 bondlength 
and would affect the shape of the curve. Using a Morse 
potential fit, we obtained an estimated bondlength of 
1.416(4) Bohr. (Using a 2nd or 4th order polynomial fit 
leads to similar results; with 4th order fit the error bar 
is three times larger). For comparison, the exact equi- 
librium bondlength of H2 is 1.40083 a.u [18], and the 
DFT/GGA bondlength is 1.4213 Bohr. 

The energy difference between the ortho- and para- 
H2 spin states (-^E and '^E, respectively) was also cal- 
culated. We note that for the singlet H2 two-electron 
system, a HS transformation based on the magnetization 
[26] can be made to eliminate the sign problem and thus 
the need for the phaseless approximation. In this case 
the AF QMC calculations will become exact. This is not 
done here, since our goal is to benchmark the general al- 
gorithm. The calculations for both ortho- and para-H2 
were at the experimental bondlength of ortho II2. Ta- 
ble II summarizes the results. The exact value obtained 
by Kolos and Roothaan is 10.495 eV [18], with which the 
AF QMC value at the larger supercell size is in excellent 
agreement. 



B. H2+H — > H+H2 transition state 

The problem of calculating the transition state of II3 
is well benchmarked using a variety of methods [27-31]. 
The activation energy for the reaction II2-I-II — > H-I-II2 
is defined as the difference between the energy of the II3 
saddle point and that of the well separated H atom and 
H2 molecule. 
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TABLE III: Symmetric coUinear H3 transition state energies 
using planewaves with pseudopotentials. Results are shown 
from density functional GGA [with (PSP) and without (AE) 
pseudopotentials, respectively], DMC, and the present AF 
QMC methods. (The "all-electron" GGA(AE) results are 
from well converged large Gaussian basis set calculations.) 
The calculated results are for the linear H3 molecule with 
Ri = R2 = R, for three values of R. All energies are in 
eV. Statistical errors are on the last digit and are shown in 
parenthesis. 



R 


GGA(AE) 


GGA(PSP) 


DMC (exact) 


AF QMC 


1.600 


0.297 


0.30 


0.543 09(8) 


0.54(3) 


1.757 


0.156 


0.16 


0.416 64(4) 


0.43(3) 


1.900 


0.222 


0.22 


0.494 39(8) 


0.48(4) 



Density functional methods are not very accurate in 
calculating the activation energy. For example, DFT 
with an LDA functional gives H3 as a bound molecule 
with a binding energy of 0.087 eV at the symmetric con- 
figuration with = R2 = 1.795 Bohr. DFT/GGA, on 
the other hand, gives a barrier height of 0.152 eV at the 
symmetric configuration R = 1.767 Bohr [31]. The ex- 
perimental barrier height is 9.7 kcal/mol= 0.420632 eV 
[32]. 

Using the AF QMC method, we studied the coUincar 
H3 system for three configurations with Ri=R2=1.600, 
1.757, and 1.900 Bohr. Table III shows the calcu- 
lated barrier heights and compares these to results from 
DFT/GGA all-electron and pseudopotential calculations, 
and to results from recent DMC calculations [29]. The 
AE and pseudopotential DFT/GGA results are in excel- 
lent agreement with each other, a further indication of 
the good quality of the H pseudopotential. The DMC 
calculations [29] are exact in this case, through the use 
of a cancellation scheme [9, 33], which is very effective 
at eliminating the sign problem for small systems. The 
AF QMC values are in good agreement with the exact 
calculated results. 

Planewave based AF QMC calculations of the H3 tran- 
sition state are very expensive, since the energy variations 
in the Born-Oppenheimer curve are quite small as seen in 
Table III. To achieve the necessary accuracy, large super- 
cells are needed, which results in large planewave basis 
sets. The large basis sets lead to many thousands of AF's 
in Eq. (3). Moreover, a large number of AF's in general 
lead to a more severe phase problem and thus potentially 
a more pronounced role for the phaseless approximation. 
The larger AF QMC statistical errors, compared to the 
highly optimized DMC results as well as to our Gaus- 
sian basis results in Sec. IV B, reflect the inefficiency of 
planewave basis sets for isolated molecules. These calcu- 
lations are valuable, despite their computational cost, as 
they demonstrate the robustness of the method. 



C. Hydrogen-bonded molecules 

Complementing the study above of the H3 system, 
where energy differences arc small, we also examined 
three other hydrogen-bonded molecules: H2O, OH, and 
H2O2, where the energy scales are large. Table IV com- 
pares the binding energies calculated using DFT/GGA 
(both pseudopotential and all-electron), DMC [34], and 
the present AF QMC method. (Results for the O2 and 
O3 molecules are included, because they are pertinent 
to the discussion of pseudopotentials errors below.) The 
experimental values [35], with the zero point energy re- 
moved, are also shown. All of the calculations are per- 
formed at the experimental geometries of the molecules. 
The density functional all-electron binding energies in 
Table IV were obtained using the highly converged triple- 
zeta AND basis sets of Widmark, Malmqvist, and Roos 
[36]. They are in good agreement with published all- 
electron results. For example, the all-electron binding 
energy of H2O is 10.147 eV and that of OH is 4.77 eV 
in Ref. [24]. In Ref. [25], the binding energy of H2O is 
10.265 eV, and that of O2 is 6.298 eV [25]. 

In all of the molecules except H2O, the DFT pseu- 
dopotential result seems to be in better agreement with 
the experimental value than the all-electron result. This 
is fortuitous and by no means suggest that the pseu- 
dopotential results are better than the all-electron values, 
since the pseudopotentials results should reproduce the 
all-electron value obtained with the same theory. Any 
differences are in fact due to transferability errors of the 
pseudopotentials. At the density functional level, the 
molecular systems H2O, OH, and H2O2 all need a non- 
linear core correction (NLCC). The NLCC was intro- 
duced into DFT pseudopotential calculations by Louie 
et. al. [22]. It arises from the DFT-generated pseudopo- 
tential for oxygen, at the pseudopotential construction 
level in the descreening step, where the valence Hartree 
and nonlinear exchange-correlation terms are subtracted 
to obtain the ionic pseudopotential. The Hartree term 
is linear in the valence charge and can be subtracted ex- 
actly. This is not the case with the nonlinear exchange- 
correlation potential, and will lead to errors especially 
when there is an overlap between the core and the va- 
lence charge densities. According to the NLCC correc- 
tion scheme, this error can be largely rectified by retain- 
ing an approximate pseudo-core charge density, and car- 
rying it properly in the target (molecular or solid) cal- 
culations. This generally improves the transferability of 
the pseudopotentials [22, 23]. The problem of NLCC is 
absent in effective core-potentials (ECP) generated using 
the Hartree-Fock method. 

All of the molecules in Table IV suffer from the 
NLCC error which originates predominantly from the 
spin-polarized oxygen atom, where the NLCC can be as 
large as 0.3 eV/atom within a GGA-PBE calculation [23]. 
For this reason, we have also included results for the O2 
and O3 molecules. (The AF QMC value for O2 is taken 
Ref. [13].) As seen in the table, the binding energies 
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TABLE IV: Calculated binding energies of H2O, OH, H2O2, 
O2, and O3. Results are shown from density functional 
GGA [with (PSP) and without (AE) pseudopotentials, re- 
spectively], DMC, and the present AF QMC methods. Ex- 
perimental results are also shown. DFT/GGA(PSP) and the 
present AF QMC results were calculated using plancwave ba- 
sis sets with pseudopotentials. DFT/GGA(AE) is calculated 
using highly converged Gaussian basis sets. The DMC [34] 
results were also obtained using pseudopotentials. The zero 
point energy is removed from the experimental data [35]. All 
energies are in eV. Statistical errors are on the last digit and 
are shown in parenthesis. 





GGA(AE) 


GGA(PSP) 


DMC 


AF QMC 


Expt. 


H20 


10.19 


9.82 


10.10(8) 


9.9(1) 


10.09 


OH 


4.79 


4.60 


4.6(1) 


4.7(1) 


4.63 


H2O2 


12.26 


11.66 


11.4(1) 


11.9(3) 


11.65 


O2 


6.22 


5.72 




5.2(1) 


5.21 


O3 


7.99 


7.12 




6.2(2) 


5.82 



of H2O, OH, H2O2, O2, and O3 are smaller than the 
corresponding all-electron values by ~ 0.37, 0.19, 0.60, 
0.50, and 0.87 cV, respectively. These values are approx- 
imately proportional to the number of oxygen atoms in 
the corresponding molecule with a proportionality con- 
stant « 0.3 eV, which agrees with the value reported in 
Ref. [23]. 

The pseudopotcntial is of course used differently in 
many-body AF QMC calculations. Despite the need 
for NLCC at the DPT level, the oxygen pseudopotcn- 
tial seems to be of good quality when used in AF QMC. 
In all the cases, the AF QMC results are in good agree- 
ment with DMC and with the experimental values. The 
largest discrepancy with experiment is « 0.4(2) cV with 
O3, and it is in opposite direction to the NLCC as done 
at the density functional level. 

The need for the non-linear core correction does not 
indicate a failure of the frozen-core approximation, but 
rather is a consequence of the non-linear dependence of 
the spin-dependent exchange-correlation potential on the 
total spin-density (valenceH-core) in density functional 
theory. The QMC calculations depend only on the bare 
ionic pseudopotcntial and do not have this explicit de- 
pendence on the (frozen) core-electron spin densities. It 
is thus reasonable to expect the QMC results to be not 
as sensitive to this issue. 



IV. RESULTS USING GAUSSIAN BASIS SETS 

In this section, we present our studies using Gaussian 
basis sets. For comparison, some of the systems are re- 
peated from the planewave and pseudopotcntial studies 
in the previous section. Gaussian basis sets arc in general 
more efficient for isolated molecules. For example, the 
calculations below on the Van der Waals minimum in H3 
would be very difficult with the plancwave formalism, be- 
cause of the large supercells necessary, and because of the 



TABLE V: Symmetric collinear H3 transition state total en- 
ergies using aug-cc-pVDZ and aug-cc-pVTZ Gaussian basis 
sets [6]. We examined 5 configurations with Ri — R2 = R, 
and we report the unrestricted Hartrcc-Fock (UHF), full con- 
figuration interaction (FCI), and AF QMC total energies. 
Bondlengths are in Bohrs and energies are in Hartrees. Statis- 
tical errors are on the last digit and are shown in parenthesis. 



R UHF 


FCI 


AF QMC 


aug-cc-pVDZ 






1.600 -1.595 026 


-1.642 820 


-1.642 56(5) 


1.700 -1.600 252 


-1.648 186 


-1.647 75(5) 


1.757 -1.601336 


-1.649 328 


-1.648 82(5) 


1.800 -1.601406 


-1.649 433 


-1.648 98(6) 


1.900 -1.599 536 


-1.647606 


-1.646 97(6) 


aug-cc-pVTZ 






1.600 -1.599 843 


-1.652 219 


-1.651 78(7) 


1.700 -1.604 162 


-1.656 405 


-1.655 86(7) 


1.757 -1.604 835 


-1.657013 


-1.656 52(7) 


1.800 -1.604 638 


-1.656 770 


-1.656 24(8) 


1.900 -1.602 269 


-1.654 285 


-1.653 68(9) 



high statistical accuracy required to distinguish the small 
energy scales. Also, all-electron calculations are feasible 
with a Gaussian basis, at least for lighter elements, so 
systematic errors due to the use of pseudopotentials can 
be avoided without incurring much additional cost. 

Direct comparison with experimental results requires 
large, well-converged basis sets in the AF QMC calcula- 
tions [3, 14]. As mentioned, the convergence of Gaussian 
basis sets is not as straightforward to control as that of 
planewaves. For benchmarking the accuracy of the AF 
QMC method, however, we can also compare with other 
established correlated methods such as full configuration 
interaction (FCI) and CCSD(T), since all the methods 
operate on the same Hilbert space. FCI energies are the 
exact results for the Hilbert space thus defined. The FCI 
method has an exponential scaling with the number of 
particles and basis size, so it is only used with small sys- 
tems. In this section, we study H2 and H3, which are 
challenging examples for mean-field methods, and com- 
pare the AF QMC results with exact results. 

The matrix elements which enter in the definition of 
the Hamiltonian of the system of Eq. (1) are calculated 
using NWCHEM [3, 17]. The trial wave functions, which 
are used to control the phase problem, are mostly com- 
puted using unrestricted Hartrce-Fock (UHF) methods, 
although we have also tested ones from density functional 
methods. In previous studies, we have rarely seen any dif- 
ference in the AF QMC results between these two types of 
trial wave functions. This is the case for most of the sys- 
tems in the present work, and only one set of results are 
reported. In H3 near the Van dcr Waals minimum, where 
extremely small energy scales need to be resolved, we find 
small differences (~ 0.1 milli-Hartree), and we report re- 
sults from the separate trial wave functions. The FCI 
calculations were performed using MOLPRO [37, 38]. 
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TABLE VI: H3 total energies in the van der Walls limit. _Ri 
is fixed at 1.4 Bohr, and R2 is varied between 4 and 10 Bohr. 
The aug-cc-pVTZ basis set is used. Energies are in Hartrees. 
Statistical errors are on the last digit and are shown in paren- 
thesis. 



R2 



4 
5 
6 
7 
10 



FCI 



-1.671 577 
-1.672 455 
-1.672 535 
-1.672 508 
-1.672 462 



AF QMC/UHF 



-1.671 60(9) 
-1.672 50(8) 
-1.672 63(6) 
-1.672 65(5) 
-1.672 57(6) 



A. Bondlength of H2 




Dissociation limit - 

14 15 



We first study H2 again, witli a cc-pVTZ basis set 
which has 28 basis functions for the molecule. This is 
to be compared with the planewave calculations which 
has about 5,000 to 70,000 planewaves for the different 
supercells used. These H-bonded systems are especially 
favorable for localized basis sets. The AF QMC equilib- 
rium bondlength R = 1.4025(6) Bohr compares well to 
the corresponding FCI bondlength of i? = 1.40265 Bohr, 
with both methods using the cc-pVTZ basis. This is 
a substantially better estimate of the exact infinite ba- 
sis result of Re = 1.40083 Bohr [18] than was obtained 
from the planewave AF QMC results in Fig. 1. The re- 
maining finite-basis error is much smaller than the sta- 
tistical errors in the planewave calculations. (The small 
residual finite-basis error is mostly removed at the cc- 
pVQZ basis set level, with an equilibrium bondlength of 
R = 1.40111 Bohr from FCL) 



FIG. 2: Potential energy curve of H3 in the van der Waals 
limit using aug-cc-pVDZ basis set. Ri is set to 1.4 Bohr, 
and R2 is varied between 4 and 10 Bohrs. (The dissociation 
limit is shown at R2 = 15 Bohr in the figure). FCI results 
are compared with AF QMC results with three different trial 
wave functions, from UHF and DFT with CCA and B3LYP 
functionals, respectively. The inset shows the correspond- 
ing potential energy curves obtained from UHF, GGA, and 
B3LYP. (For clarity, the UHF and GGA energies are shifted 
by —0.047 and —0.154 Hartrees in the inset, respectively.) 



only a small fraction of the computational time. Even 
with these relatively small basis sets, we see that the 
finite-basis errors are quite small here. In fact, the FCI 
barrier height with the aug-cc-pVTZ basis is in agree- 
ment with the experimental value 0.420632 eV [32]. 



B. H2+H — > H+H2 transition state 

Table V presents calculated total energies of H3 with 
aug-cc-pVDZ and aug-cc-pVTZ basis sets [6]. Results ob- 
tained with UHF, FCI, and the present AF QMC meth- 
ods are shown for five different geometries in the coUinear 
H3 system. The present TZ-basis FCI results were cross- 
checked with those in Ref . [30] , which contains a detailed 
study of the Born-Oppenheimer potential energy curves 
for the H-I-H2 system. 

The AF QMC total energies are in excellent agree- 
ment, to within less than 1 uiEh, with the FCI energies. 
The AF QMC barrier heights with the aug-cc-pVDZ and 
aug-cc-pVTZ basis sets at i? = 1.757 Bohr are 0.444(2) 
and 0.434(3) eV, respectively. The corresponding FCI 
results are 0.4309 and 0.4202 eV, respectively. Thus the 
AF QMC results show a systematic error of ^ 0.015 eV 
in the barrier height. It is possible to resolve these small 
discrepancies, because the basis sets are much more com- 
pact, with 25 to 75 Gaussian basis functions as opposed 
to approximately 10, 000 planewaves in the calculations 
in Sec. IIIB. As a result, the statistical errors are smaller 
than in the planewave calculations by a factor of 10, with 



C. Van der Waals minimum in collinear H3 

The van der Waals minimum of II3 is studied by fixing 
i?i = 1.4 Bohr (the H2 equilibrium bond length), while 
the distance R2 between the third H atom and the closer 
of the two atoms in II2 was varied between 4 and 10 
Bohrs. The potential energy curve of this system exhibits 
a very shallow minimum of approximately 85 iiEh [30] at 
i?2 ~ 6 Bohr. Two different basis sets, aug-cc-pVDZ and 
aug-cc-pVTZ, were used. Table VI shows the aug-cc- 
pVTZ results, and Fig. 2 plots the aug-cc-pVDZ results. 

As seen in Table VI, the AF QMC all-electron total en- 
ergies are in excellent agreement with FCI, with a maxi- 
mum discrepancy of about 0.14(5) mEjy. The AF QMC 
energies, which are calculated with the mixed-estimator, 
are not variational, as is evident in the results from both 
basis sets compared to FCI. The AF QMC results in Ta- 
ble VI are obtained with an UHF trial wave function. In 
most of our molecular calculations with Gaussian basis 
sets, the UHF solution, which is the variationally opti- 
mal single Slater determinant, has been chosen as the 
trial wave function [3, 14]. In the present case, the UHF 
method actually fails to give a Van der Waals minimum, 
as can be seen from the inset of Fig. 2. It is reassuring 
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that AF QMC correctly reproduces the minimum with 
UHF as a trial wave function. 

The effects of using two other single Slater determi- 
nant trial wave functions were also tested. These were 
obtained from DFT GGA and B3LYP calculations with 
the aug-cc-pVDZ basis set. The corresponding results are 
also shown in Fig. 2. In the DFT calculations (shown in 
the inset of Fig. 2), both GGA and B3LYP predict the ex- 
istence of a minimum, although B3LYP gives an unphys- 
ical small barrier at about R2 « 7 Bohr. The AF QMC 
results obtained with UHF, GGA, and B3LYP Slater de- 
terminants as trial wave functions differ somewhat, but 
arc reasonably close to each another. With the GGA trial 
wave function, AF QMC "repairs" the well depth (possi- 
bly with a slight over-correction). With the B3LYP trial 
wave function, AF QMC appears to under-estimate the 
well-depth, giving a well-shape that is difficult to charac- 
terize because of the statistical errors and the extremely 
small energy scale of these features. 

V. SUMMARY 

We have presented a benchmark study of the phaseless 
AF QMC method in various H-bonded molecules. The 
auxiliary-field QMC method is a many-body approach 
formulated in a Hilbcrt space defined by a single-particle 
basis. The choice of a basis set is often of key impor- 
tance, as it can affect the efficiency of the calculation. 
In the case of AF QMC, the basis set choice can also 
affect the systematic error, because of the different HS 
transformation that can result. In this study, we em- 
ployed planewave basis sets with pseudopotentials and 
all-electron Gaussian basis sets, to compare the perfor- 
mance of the AF QMC method. The planewave HS de- 
composition was tailored to the planewave representa- 
tion, resulting in 0(8 M) auxiliary fields, where M is 
the number of planewaves. For the Gaussian basis sets, 
the generic HS decomposition described in Section II was 
used, resulting in C(M^) atixiliary fields. Typical M val- 
ues in this study were tens of thousands in the planewave 
calculations and a hundred in the Gaussian calculations. 

The planewave calculations were carried out for H2, 
H2 -|- H near the transition state, H2O, OH, and H202. 
Non-linear core corrections to the oxygen pseudopoten- 
tial were discussed using additional calculations for the 
O2 and O3 molecules. DFT GGA pseudopotentials were 
employed. The trial wave functions were single Slater 
determinants obtained from DFT GGA with identical 
planewave and pseudopotential parameters as in the AF 
QMC calculations. Hard pseudopotentials and large 



planewave cutoffs were used to ensure basis-size con- 
vergence and the transferability of the pseudopotentials. 
Large supercells were employed to remove finite-size er- 
rors. To mimic typical systems in the solid state, no op- 
timization was done to take advantage of the simplicity 
of these particular systems. The binding energies com- 
puted from AF QMC have statistical errors of 0.1-0.3 eV 
as a result. Within this accuracy, the AF QMC results 
are in excellent agreement with experimental values. 

Gaussian basis AF QMC calculations were carried out 
on H2, the transition state of H2 -|- H, as well as the 
van der Waals minimum in linear H2 + H. These calcu- 
lations are within the framework of standard quantum 
chemistry many-body using the full Hamiltonian with- 
out pseudopotentials. UHF single Slater determinants 
were used as the trial wave function. For various geome- 
tries, the absolute total energies from AF QMC agree 
with FCI to well within 1 ibEh- The calculated equilib- 
rium bondlengths and potential energy curves are also in 
excellent agreement with FCI. In H2 -I- H, AF QMC cor- 
rectly recovers the van der Waals well with a UHF trial 
wave function which in itself predicts no binding. 

Comparing planewave and Gaussian basis set AF QMC 
results, we can conclude the following. In the Gaussian 
basis calculations, as evident from FCI comparisons, er- 
rors due to controlling the phase problem in the phaseless 
approximation are well within 1 irEh in the absolute en- 
ergies. Achieving the infinite basis limit is more straight- 
forward using planewave based AF QMC, but statistical 
errors are larger for the isolated molecules studied due 
to the need for large supercells. Within statistical errors, 
however, the AF QMC results using both types of basis 
sets were in agreement. This indicates that errors due 
to the use of pseudopotentials with planewave basis sets 
were smaller than the statistical errors. Finally, within 
statistical errors, the performance of the phaseless AF 
QMC method, did not appear to be sensitive to the type 
of HS decomposition used, despite drastic differences in 
basis size and the number of auxiliary fields, was tailored 
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